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Abstract 

We introduce an adaptation algorithm by which an ensemble of coupled oscillators 
with attractive and repulsive interactions is induced to adopt a prescribed syn- 
chronized state. While the performance of adaptation is controlled by measuring a 
macroscopic quantity, which characterizes the achieved degree of synchronization, 
adaptive changes are introduced at the microscopic level of the interaction network, 
by modifying the configuration of repulsive interactions. This scheme emulates the 
distinct levels of selection and mutation in biological evolution and learning. 
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1 Introduction 

The emergence of coherent collective behaviour, out of the interaction of a 
large number of active elements with relatively simple individual dynamics, is 
the key signature of a complex system [T] . The nature of the mechanisms which 
underly such emergence, however, can substantially vary between different sys- 
tems. In physical or chemical phenomena -for instance, in the formation of 
spatiotemporal patterns in reacting fluids- collective behaviour is the unavoid- 
able, spontaneous consequence of elementary mechanical processes occurring 
at microscopic scales. In biological systems, on the other hand, collective be- 
haviour is not just the macroscopic manifestation of the mutual organization of 
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microscopic entities. It is, as well, the consequence of a much longer evolution- 
ary process, by which certain macroscopic traits are selected at the expense of 
others, while they change from one generation to the next according the laws 
of heredity. The joint action of natural selection and mutations, as described 
by the Darwinian theory of evolution, makes it possible that, in a given envi- 
ronment, the collective performance of a class of biological complex systems 
-a species- improves as generations succeed each other. Similarly, the acquisi- 
tion of certain forms of behaviour by learning implies a process by which the 
response of a living being to given stimuli becomes specific and optimized. 

It is interesting that -both in evolution and in learning- the mechanism of 
selection, on one hand, and the changes that eventually improve performance, 
on the other, act at very distinct levels. Natural selection of beneficial traits, 
as well as evaluation of learnt tasks, take place at macroscopic scales, at the 
level of the organic interaction of individuals with their environment. Variation 
processes, in contrast, take place at microscopic scales, genetic for evolution 
and neural for learning. It is precisely the emergent nature of macroscopic 
traits which connects the two levels. 

In this paper, we explore the possibility of emulating this scenario with a rather 
simple model, based on the collective dynamics of coupled oscillators. It is well 
known that ensembles of interacting oscillators may undergo self-organization 
processes which lead to different forms of synchronized behaviour [2^3j. A vari- 
ety of collective states -in the form of single or multiple synchronized clusters, 
for instance- can be realized by altering the pattern of interactions between 
oscillators. We exploit this behavioral diversity for inducing an oscillator en- 
semble to adopt a specific, target synchronization state, through adaptive 
introduction of gradual changes in its interaction pattern. Imitating biological 
adaptation, achievement of the goal is evaluated through a macroscopic mea- 
sure of synchronization, while changes are introduced at the microscopic level 
of pair interactions. 

We analyze, in the next section, the possible synchronization states of an oscil- 
lator ensemble with both attractive and repulsive interactions. We show that 
if repulsive interactions are conveniently distributed, the ensemble splits into 
two synchronized clusters with opposite phases. Choosing this two-cluster con- 
figuration as the target state, in Section [3] we apply an adaptation algorithm 
to induce an ensemble, whose repulsive interactions are initially distributed 
at random, to modify its interaction pattern in such a way that its collec- 
tive behaviour approaches the target. This formulation contrasts with previ- 
ous studies of network adaptation to self-organization |3], where the degree of 
synchronization is evaluated at local -not collective- level. In spite of the huge 
number of possible ways of distributing the repulsive interactions, our system 
manages to satisfactorily achieve its goal. 
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2 Oscillator networks with repulsive interactions 



We consider a Kuramoto-like model for coupled phase oscillators, governed 
by the equations 



K ^ 



where G [0, 27r) is the phase of oscillator i, and Ui is its natural fre- 

quency. The positive constant K measures the coupling between oscillators, 
and Wij = ibl defines the sign of the interaction between oscillator i and 
j. Positive and negative Wij represent, respectively, attractive and repulsive 
"forces" between oscillators. Interactions are assumed to be symmetric, so 
that Wij = Wji. In the following, it will be useful to conceive the interaction 
pattern as defining a network with vertices occupied by the oscillators, whose 
links join oscillators with attractive interactions. The complementary network, 
with links corresponding to repulsive interactions, provides an alternative rep- 
resentation [6]. 

It is well-known that in the case of global attractive coupling, Wij = 1 for all 
i and j, and in the thermodynamical limit, N ^ oo, the oscillator ensemble 
undergoes a synchronization transition as the coupling strength K grows [2]. 
This transition is quantified by the Kuramoto order parameter 




where (■) indicates averages over sufficiently long times. Below the critical 
coupling strength Kc, phases are uniformly distributed in [0, 27r), so that z = 0. 
Above Kc, on the other hand, z grows monotonically with K, revealing a 
progressive condensation of phases. At the same time, an increasingly large 
cluster of synchronized oscillators, whose time-averaged frequencies, 

1 T 

rii = {<P,) = lirn^- J ^i{t)dt, (3) 





are mutually identical, dominates the collective dynamics of the ensemble. 
The critical point Kc depends on the distribution of natural frequencies cUj: 
a higher dispersion in the natural frequencies requires a stronger coupling to 
synchronize the ensemble. For asymptotically large K, the cluster of synchro- 
nized oscillators entrains the whole ensemble, and all the phases collapse to 
the same value, so that z ^ 1. 
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As it may be expected, when some of the interactions are repulsive, synchro- 
nization is harder to attain. This was originally confirmed by Daido [7j, who 
studied Eqs. ([T]) in the case where the symmetric weights Wij are drawn at 
random from a Gaussian distribution centered at zero. Moreover, it turns out 
that the emerging patterns of collective behavior depend crucially on the form 
in which repulsive interactions are distributed over the ensemble. In Ref. [6], 
the effect of a random (Erdos-Renyi 0) network of repulsive interactions was 
studied for identical phase oscillators (cUj = a; for all i). It was shown that, 
as the number of repulsive interactions grows, fully synchronized ensembles 
undergo a sharp transition to desynchronization. Figure [U on the other hand, 
shows results for ensembles of = 10^ nonidentical oscillaXors with an Erdos- 
Renyi network of repulsive interactions. Natural frequencies are drawn at ran- 
dom from a Gaussian distribution g{uj) = exp(— co'^/2)/-\/27r. Each interaction 
weight is given by 



so that X is the average fraction of repulsive interactions. The plot displays the 
order parameter z, Eq. ([2]), as a function of the coupling strength K for several 
values of the fraction x. Note that a; < 1/2 in all cases, so that the number of 
attractive interactions is always larger than that of repulsive interactions. 

From Kuramoto's theory p], it is known that in the absence of repulsive inter- 
actions (x = 0) the synchronization transition for the present distribution of 
natural frequencies takes place at Kc ~ 1.6. Our results show that, as x grows, 
the transition shifts to higher values of K, and the transition zone becomes 
broader. Nevertheless, for sufficiently large coupling intensities, the order pa- 
rameter seems to approach its maximum value, z = 1, which corresponds 
to the state of full synchronization. Thus, as coupling becomes stronger, the 
more abundant attractive interactions overcome the effect of repulsive interac- 
tions, and the whole ensemble aggregates into a single cluster with well-defined 
phase. 

Other modes of collective behaviour can be expected if, instead of being ran- 
dom, the networks of attractive and repulsive interactions are given some spe- 
cial structure. Specifically, consider that the oscillator ensemble is divided into 
two groups, such that inside each group all interactions are attractive, whereas 
they are repulsive between oscillators of different groups. In this situation, the 
network of attractive interactions is fully connected inside each group, and the 




1 with probability 1 — x 



1 with probability x. 



(4) 
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Fig. 1. The Kuramoto order parameter 2; as a function of the coupUng strength 
K, for an ensemble of 10"^ oscillators with a fraction x of repulsive interactions 
distributed at random. Repulsive interactions define an Erdos-Renyi network over 
the oscillator ensemble. 

two groups are mutually disconnected. The interaction weights are given by 

{1, if i and j belong to the same group, 
(5) 
— 1, otherwise. 



If the groups have, respectively, A^i and N2 oscillators (A^i + A^2 = N), the 
fraction of repulsive interactions is x = 2N1N2/ N{N — 1), or equivalently. 



2" 



1 ± ^l-2x{N -l)/N ^ — (1 ± Vl - 2x) 



(6) 



where the approximation holds for large A^. Note that these numbers make 
sense if x < N/2{N - 1) ^ 1/2. 

It can be easily realized that, if all natural frequencies are equal {uJi = uj for all 
i), the only stable solution to Eqs. ([T]) with the interactions given by Eq. ^ 
corresponds to all the oscillators of each group having exactly the same phase, 
and the two groups having opposite phases. In other words, an ensemble of 
identical oscillators with the interactions of Eq. ([5]) splits into two point-like 
opposite clusters in the phase coordinate, each cluster corresponding to one 



group. When natural frequencies are not mutually identical, we expect that 
the clusters spread in phase but, for large coupling strengths, the division into 
two opposite phase clusters corresponding to the groups should persist. This 
is in fact what numerical results show. Figure [2] displays results for the order 
parameter function of the coupling strength K for several values of x, 

with the distribution of interactions given by Eq. ([5]), in a system of = 10^ 
oscillators. Note that, now, the transition to synchronized behaviour occurs 
always at the same value of K. On the other hand, the asymptotic value of z 
for large K depends on x. Indeed, as discussed above, for large K the ensemble 
is organized in two point-like clusters of opposite phase so that, according to 
Eq. ([2]), the order parameter equals 



z{K — > oo) 



IAWV2 



2x(iV - \)IN ^ Vl - 2x, 



(7) 



an analytical prediction which is in full agreement with numerical results. 
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Fig. 2. The Kuramoto order parameter z as a function of the coupling strength 
for an ensemble of 10^ oscillators divided into two groups. The interaction between 
two oscillators is attractive (repulsive) if they (do not) belong to the same group. 
The sizes A^i and N2 of the two groups are given, as a function of the fraction x of 
repulsive interactions, by Eq. ([6]). 
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The splitting of the ensemble into two clusters is quantitatively revealed by 
an additional order parameter |9], 




[cf. Eq. (|2])], which vanishes if oscillator phases are uniformly distributed in 
[0,27r), and attains its maximum, Z2 = 1, if oscillators form two point-like 
clusters of opposite phases. In this limit, Z2 is independent of the relative size 
of the two clusters, so that Z2 = 1 also in the extreme case of a single point 
cluster -equivalent to A^^i = and N2 = 0, or vice versa. Therefore, for a 
random interaction network, both z and Z2 approach unity as the coupling 
strength grows. It can be easily shown that, in this case, z > Z2- For the 
two-group network, on the other hand, Z2 tends to unity, while z approaches 
the asymptotic value quoted in Eq. ([7]) so that, for large K, z < Z2. This is 
illustrated in Fig. [31 for a fixed value of the fraction x of repulsive interactions. 
The combined analysis of z and Z2 as a function of K thus discerns between the 
two different forms of collective behaviour of the oscillator ensemble, derived 
from the different structures of their interaction networks. 
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Fig. 3. The order parameters z and Z2 as functions of the coupling strength K, 
for an ensemble of 50 oscillators with a fraction x = 0.3 of repulsive interactions 
distributed at random (dashed lines) or between two groups (full lines). 
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3 Adaptation to collective organization 



In the preceding section, we have shown that a specifically structured interac- 
tion network induces splitting of the oscillator ensemble into two synchronized 
clusters with opposite phases. In contrast, a random distribution of attractive 
and repulsive interactions gives rise to a single synchronized cluster. Suppose 
now to have an ensemble with randomly distributed interactions. Would it 
be possible, following our discussion in the Introduction, to adaptively evolve 
the interaction pattern toward the two-group structure by applying selection 
pressure at the level of collective behaviour? This would amount to modifying 
microscopic attributes by selection of macroscopic features, as in biological 
evolution and learning. 

In order to explore this possibility, we implement an algorithm of adaptation of 
the interaction pattern based on a Monte Carlo scheme. The algorithm is con- 
trolled by monitoring the collective dynamical state of the oscillator ensemble, 
as quantified by the two-cluster order parameter Z2- We work at intermediate 
values of the coupling strength, where the difference in Z2 between the two dy- 
namical modes studied above is larger (cf. Fig. [3]). In this intermediate region, 
coupling is sufficiently strong to induce synchronized collective behaviour, but 
does not reach the regime where oscillators are confined to almost point-like 
clusters. In our calculations, we choose K = A. 

We start from an interaction pattern with a fraction x of repulsive interactions 
distributed at random all over the ensemble, with weights Wij. Equations ([1]) 
are solved by means of an Euler algorithm, with time increment h = 0.1. The 
system is left to evolve for a few hundred time units (typically, T = 250) so 
that Z2 attains a well-defined value. Then, a small change is introduced in the 
interaction pattern. Two oscillator pairs, {k,l} and {r,s}, with interactions 
of opposite signs, W^i = 1 and Wrs = —1, are chosen at random and their 
interaction weights are interchanged, so that Wj^i = — 1 and W^g = 1. This 
procedure is repeated on np randomly chosen oscillator pairs. With this new 
interaction pattern, the system evolves again and a new value of the two- 
cluster order parameter, Z2, is determined. If Az2 = Z2 — Z2 > 0, the change 
in the interaction weights is accepted. If on the other hand, Az2 < 0, the 
change is accepted with probability exp(Az2/^), and rejected otherwise. The 
fictitious temperature 9 is used to control the speed of the adaptation process. 
The algorithm is repeated a number of times (typically, several hundred), 
until the successive values of Z2 converge to a stationary level. Following the 
standard Monte Carlo practice, the temperature 9 is slowly decreased as the 
process progresses, to gradually attenuate the effect of fluctuations in Z2- In 
our simulations, the number of pairs whose interaction weights are modified 
at each adaptation step is around np = 10, and is adjusted to make the final 
value of Z2 as large as possible. 
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To illustrate our results, we consider a realization of the above Monte Carlo 
process in an ensemble of iV = 50 oscillators, where the fraction of repulsive 
interactions is x = 3/7 ~ 0.43. In a two-group interaction network, this would 
correspond to groups of sizes A^^i = 15 and N2 = 35. From an initial value 
Z2 ~ 0.14, obtained for a random interaction network, the two-cluster order 
parameter significantly increases by a factor slightly above four, reaching Z2 ~ 
0.57. Numerical results for systems of this size, with K = 4 and the two-group 
structure given by the interaction weights of Eq. ([5]), show however that the 
two-cluster order parameter can reach values around Z2 = 0.9. After the Monte 
Carlo process, the obtained value of Z2 is still considerably below this level. In 
any case, in order to evaluate to which extent has the system evolved towards 
the target state, we examine both its resulting dynamical behaviour and its 
interaction pattern. 




Fig. 4. Number of oscillator pairs {i,j} with relative phase Acpij = |0j ~ and 
attractive or repulsive interactions. Individual phases were measured at a fixed time 
after adaptation to the two-group state, in a realization with the parameters quoted 
in the text. 



3.1 Distribution of phases and synchronization 



Figure m shows results illustrating the correlation between the relative phase of 
each oscillator pair and the sign of the corresponding interaction, after adap- 
tation. For two oscillators i and j, the relative phase A0jj = |0i — 0j| should 
be close to or tt for Wij = 1 or — 1, respectively. Indeed, as expected, the his- 
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tograms of Fig. H] display a clear excess of attractive (repulsive) interactions 
for small (large) relative phases. A quantitative measure of the correlation 
between A0jj and Wij is given by the index 

C= ^^J^_-^^ I:W^,cosA<|>,,. (9) 



It is easily shown that, for the synchronized state with two point clusters - 
which would be achieved for large coupling strengths and a perfect two-group 
interaction network- we get C = 1. When oscillators are uniformly spread in 
phases, on the other hand, C ~ A^~^. The latter is the expected value of C for 
a random interaction network. In fact, from Fig.[T], we note that for the present 
values of x and K the oscillator ensemble with randomly distributed repulsive 
interactions is essentially unsynchronized. For the data corresponding to Fig. 
m in contrast, we find C ~ 0.22. As may have been expected, due to the 
moderate coupling strength, this is still far from the maximum. However, C is 
an order of magnitude larger than the value expected for an unsynchronized 
ensemble. This indicates that adaptation has had a substantial effect in the 
desired direction. 




Fig. 5. Two snapshots of the distribution of phases after adaptation to the two-group 
state, in an ensemble of 50 oscillators with the parameters specified in the main text. 
Individual phases are represented over the unit circle, and lines stand for repulsive 
interactions between oscillators. 

Closer inspection of the distribution of phases reveals a complementary conse- 
quence of the moderate value of K. The two plots of Fig. [S] are representations 
of the individual phases of the 50 oscillators over the unit circle. Lines joining 
oscillators represent repulsive interactions. Each plot corresponds to a differ- 
ent time, i.e. it shows a snapshot of the distribution of phases. In Fig. |5^, it is 
clear that most of the ensemble has split into two clusters with approximately 
opposite phases, with most of the repulsive interactions standing between os- 
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cillators in different clusters. Some time later, as illustrated by the snapshot 
of Fig. [5]d, clusters are however much less defined and the distribution over 
the unit circle is more uniform. While most of the repulsive interactions still 
correspond to pairs of oscillators with opposite phases, it is clear that the 
ensemble is considerably more spread. 




Fig. 6. Main plot: The average frequency fij as a function the natural frequency uJi 
for individual oscillators after adaptation, with the parameters specified in the main 
text. Inserts: Snapshots of the distribution of phases (cf. Fig. E]) for the synchronized 
clusters inside dotted frames. 

To disclose the origin of this apparent time variation in the degree of organi- 
zation in phases, it is useful to consider as well the distribution of frequencies 
over the ensemble. The main plot in Fig.[6]represents the average frequency 
defined in Eq. ([3]), as a function of the natural frequency Ui for each individual 
oscillator. In this kind of plot, horizontal arrays of dots represent clusters of 
oscillators whose natural frequencies are different, but which have synchro- 
nized to a common average frequency. At the present value of K, we find 
several of these clusters, of different sizes, with different average frequencies. 
The analysis of the distribution of repulsive interactions inside each synchro- 
nized cluster shows that the adaptation algorithm has in fact succeeded at 
separating groups with opposite phases. This is displayed for two of the syn- 
chronized clusters in the inserts of Fig. [61 However, since each synchronized 
cluster has its own average frequency, they move around the unit circle with 
respect to each other. Thus, their relative position changes, and the clustering 
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in the distribution of phases is at times more evident, at times less defined. 
This observation brings to hght a subtle difference between the goal of the 
adaptation algorithm and the synchronization of the ensemble: for intermedi- 
ate values of the coupling strength, pairs of clusters with opposite phases form 
due to the redistribution of repulsive interactions, though desynchronization 
between different pairs may persist. 



3.2 Structure of the interaction network 

The results discussed in Section 13.11 allow us to appraise the performance 
of adaptation at the level of the dynamical state achieved by the ensemble, 
through the correlation between phases and the distribution of repulsive in- 
teractions, and through the degree of synchronization. This macroscopic view- 
point can be complemented by an evaluation of the microscopic structure of 
the resulting interaction network. In other words, we aim at evaluating to 
which extent the redistribution of interactions has transformed the initial ran- 
dom structure into the two-group target. However, if such evaluation is purely 
based on an analysis of the network structure, the resulting distribution of 
phases must be disregarded. This has the drawback that we do not know a 
priori which oscillators belong to each group. 

To tackle this problem we have adapted two algorithms of community de- 
tection in networks [TO|ITT] . whose aim is to identify groups of network nodes 
which are internally best connected. In our case, this task is translated into the 
identification of two groups which maximize the number of attractive interac- 
tions inside each group. The result is given as an ordering of the oscillators, 
with lower ranks corresponding to one group and higher ranks to the other. 
In the first method, called spectral partitioning algorithm [10], the size of the 
two groups is fixed -in our case, to A^i = 15 and N2 = 35. In the second, which 
is based on the maximization of the network modularity pT], the group size 
is free. 

Figure [7^ is a graphical representation of the adjacency matrix of the network 
of attractive interactions with the numbering resulting from the first method. 
Squares represent oscillator pairs with attractive interactions. In the ideal case 
where both the adaptation algorithm and the ordering method are perfectly 
successful, all squares should be concentrated in two non-overlapping blocks 
of sizes Ni X Ni and N2 x N2, situated along the diagonal. Though our result 
is far from this ideal configuration, a larger density of squares is apparent over 
the two blocks. Figure presents the result of using the second method, with 
free group sizes. The concentration in the expected zones is again clear, though 
the resulting groups are of shghtly different sizes: A''i ^ 20 and N2 ~ 30. 
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Fig. 7. Graphical representation of the adjacency matrices of the network of attrac- 
tive interactions, after detection of the two groups, (a) by the algorithm of spectral 
partitioning and (b) by maximizing the network modularity. 



In order to quantify the performance of adaptation at the level of the network 
structure, we may use the results depicted in Fig. [7] to introduce an index 
which measures to what extent has the adjacency matrix achieved the ideal 
division into two blocks. Let / be the fraction of sites inside the blocks (ex- 
cluding the diagonals) which are occupied by squares after adaptation, and 
fji = l — x the fraction of occupied sites for a random distribution of attrac- 
tive interactions. Note that coincides, up to random fluctuations, with the 
fraction of occupied sites inside the blocks in the initial condition. We define 



D = ^ — ^. (10) 

l-fR ^ ^ 



In the ideal final state, we would have D = 1, because the filling of the 
two blocks would be perfect (/ = 1). Meanwhile, for randomly distributed 
interactions, / ^ and D ^ 0. Taking into account the effect of fluctuations, 
in fact, we find D ~ A^~^. The values of D calculated for the results shown 
in Fig. [7] are (a) D = 0.26 and (b) 0.33. As for the case of the index C, 
Eq. (Q, these values are considerably smaller than the ideal. However, during 
adaptation, they have grown by an order of magnitude, which represents an 
important improvement towards the target interaction pattern. 
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4 Discussion and conclusion 



For the 50-oscillator ensemble, we have repeated the analysis with groups of 
various sizes {Ni = 10, 20, 25) and a different coupling strength {K = 8), 
and have always obtained consistent results -not presented here for concise- 
ness. The degree of adaptation both in dynamics and in network structure, 
as measured by the quantities C and D, was similar to that presented above. 
We have also tested the adaptation algorithm using the difference ( = z^ — z^, 
instead of the two-cluster order parameter as the macroscopic measure of 
collective organization. In view of the results displayed in Fig. [3l this criterion 
is expected to exhibit better performance for large values of K, where the 
difference between the two order parameters z and Z2 is large for the two- 
group configuration, and small for the random interaction network. Finally, 
we have verified that adding a moderate number of oscillators with random 
interactions, once the original ensemble has already undergone the adaptation 
process, does not produce a substantial change in the collective dynamical 
state of the system. All these tests point at the robustness of the results pre- 
sented in Section [31 

Regarding the collective dynamical state achieved by the ensemble after the 
adaptation process, it is interesting to emphasize that, at the moderate cou- 
pling strength considered here, a difference results between the target state 
of two clusters with opposite phases, on one hand, and ensemble synchroniza- 
tion, on the other. As a result of adaptation, in fact, most oscillators end in 
pairs of clusters with opposite phases -indicating that an important fraction 
of repulsive interactions has been redistributed as expected. However, differ- 
ent pairs of clusters have different average frequencies and, therefore, are not 
mutually synchronized. In other words, there is a certain degree of indepen- 
dence between synchronized clusters and opposite-phase clusters, as depicted 
in Fig. [6l This "multilevel" clustering may represent a rich regime from the 
viewpoint of the collective dynamics [12], and deserves further study. 

As for the network structure, the results shown in Fig. [TJ along with the 
values obtained for the quantity D, may seem a modest achievement for a 
Monte Carlo adaptation algorithm on a 50-oscillator ensemble. One should 
take into account, however, that the algorithm explores the space of possible 
distributions of N1N2 repulsive interactions over the N[N — l)/2 oscillator 
pairs. For A^^i = 15 and A^2 = 35, identifying the ideal two-group configuration 
amounts at singling a particular state out of 4.7 x 10^^^ possibilities! This huge 
number approximately equals 4^*^°, which indicates that the task is comparable 
to individualizing a specific 600-nucleotide long genetic message. Note that 
this figure is not far below the length of some short coding sequences in real 
genomes [13]. 
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We have shown in this paper that the distinctive scale separation between the 
mechanisms of mutation and selection which characterize biological evolution 
-or, equivalently, neural and supervising mechanisms in learning processes- 
can be emulated by an ensemble of very simple interacting dynamical systems, 
which adapts to a prescribed form of collective behaviour by gradual changes 
in its interaction pattern. Admittedly, the two-cluster synchronization state 
here chosen as the adaptation target is by no means as complex as some of the 
most elementary functional configurations of biological entities. However, as 
long as a quantification of the adaptation level -analogous to our two-cluster 
order parameter is identified, the same algorithm could be in principle im- 
plemented for evolution towards more sophisticated collective dynamics. The 
same consideration can be applied to systems whose elementary components 
are dynamically richer that the phase oscillators studied here, such in the 
case of chaotic units or whose interaction pattern is structurally closer to 
those observed in real biological networks -typically, much more sparse than 
our fully connected graphs. Additional dynamical ingredients may include, for 
instance, the effects of internal fluctuations and noise. 

In any case, the quantitative extent and limitations of the correspondence be- 
tween biological evolution and adaptation in networks of dynamical elements 
-in particular, with respect to the performance of the adaptation algorithm- 
as well as the possibility of implementation in other complex systems with 
different forms of self-organizing collective dynamics, is worth considering in 
future work. 



References 

[1] A. S. Mikhailov and V. Calenbuhr, From Cells to Societies. Models of Complex 
Coherent Action, Springer, Berlin, 2002. 

[2] Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence, Springer, Berlin, 
1984. 

[3] S. C. Manrubia, A. S. Mikhailov, and D. H. Zanette, Emergence of Dynamical 
Order. Synchronization Phenomena in Complex Systems, World Scientific, 
Singapore, 2004. 

[4] C. Zhou and J. Kurths, Phys. Rev. Lett. 96 (2006) 164102. 

[5] H. Daido, Prog. Theor. Phys. 77 (1987) 622. 

[6] D. H. Zanette, Europhys. Lett. 72 (2005) 190. 

[7] H. Daido, Phys. Rev. Lett. 68 (1992) 1073. 

[8] P. Erdos and A. Renyi, Publ. Math. (Debrecen) 6 (1959) 290; 
R. Albert and L. Barabasi, Rev. Mod. Phys. 74 (2002) 47. 



15 



[9] H. Daido, Prog. Theor. Phys. 88 (1992) 1213. 

[10] M. Fielder, Czech. Math. J 23 (1973) 298. 

[11] M. E. J. Newman, Phys. Rev. E 74 (2006) 036104. 

[12] A. S. Mikhailov, D. H. Zanette, Y. M. Zhai, I. Z. Kiss, and J. L. Hudson, Proc. 
Nat. Acad. Sci. USA 101 (2004) 10890. 

[13] W. S. Klug and M. R. Cummings, Genetics, Prentice Hall, Upper Saddle River, 
1994. 



16 



